Intesity of m6A-Tracer around Ki67-marked nucleoli

Introduction

The goal of this analysis is to determine the enrichment of m6A-Tracer signal around nucleoli (marked by Ki67), following Ki67 pA-DamID. Dam-only signal is the control.

Setup

Load the required libraries.

# Load libraries
library(ggplot2)
library(tidyverse)
library(ggbeeswarm)
library(pixmap)

# Input data
results_dir <- c("../ts201008_E1304_E1321_confocal_various/ts201008_pADamID_m6ATracer_analysis/",
                 "../ts201102_E1349_confocal_synchronization/ts201102_RPE_m6ATracer/ts201102_RPE_timecourse_analysis/",
                 "../ts210308_E1499_confocal/ts210316_RPE_m6ATracer/ts210308_RPE_m6ATracer_analysis/")
                 ####"../ts210402_E1546_confocal_Ki67_antibodies/ts210402_HCT116_wt_pADamID_Ki67_antibodies/ts210402_HCT116_wt_pADamID_Ki67_antibodies_analysis/")

# Prepare output directory
output_dir <- "ts220503_nucleolar_enrichment_cell_types"
dir.create(output_dir, showWarnings = FALSE)

# Prepare output files

library(knitr)
opts_chunk$set(dev=c('png', 'pdf'), fig.path = file.path(output_dir, "figures/"),
               cache = T)
pdf.options(useDingbats = FALSE)

ReadPGN <- function(file.name) {
  
  # Get image, only grey values
  suppressWarnings(read.pnm(file.name)@grey*255)
  
}

Strategy

The main problem is substantial difference in microscopy slide quality between the replicates / cell types. This results in consistent results, but also in qualitative variation. I will explore a bit whether I can solve this, but otherwise I can always just “explain” this in the text. Solving means tweaking the parameters to segment nucleoli, mostly.

1) Prepare metadata

List the different samples.

metadata <- tibble(cell_statistics = dir(results_dir, recursive = T, 
                                         full.names = T,
                                         pattern = "cells_statistics")) %>%
  mutate(sample = str_remove(basename(cell_statistics),
                             "_cells_.*"),
         sample = str_remove(sample, "r._"),
         sample = str_replace(sample, "unsync", "RPE")) %>%
  filter(! grepl("Act", sample) &
           ! grepl("Preview", sample) &
           ! grepl("t\\d+h+_", sample) &
           ! grepl("NOV|HPA", sample)) %>%
  mutate(replicate = case_when(grepl("r1_", cell_statistics) ~ "r2",
                               grepl("r2_", cell_statistics) ~ "r3",
                               T ~ "r1"),
         replicate = factor(replicate, levels = c("r1", "r2", "r3")),
         cell = case_when(grepl("HCT116", sample) ~ "HCT116",
                          grepl("K562", sample) ~ "K562",
                          T ~ "RPE"),
         cell = factor(cell, levels = c("RPE", "HCT116", "K562")),
         condition = case_when(grepl("DMSO", cell_statistics) ~ "DMSO",
                               T ~ "wt"),
         target = ifelse(grepl("Ki67", sample), "Ki67", "Dam"),
         target = factor(target, levels = c("Dam", "Ki67")),
         slide = str_remove(sample, ".*_")) %>%
  mutate(sample = paste(cell, target, condition, replicate, slide, sep = "_")) %>%
  mutate(objects = str_replace(cell_statistics, 
                               "cells_statistics.csv",
                               "dapi_segment.pgm"),
         cell_internal = str_replace(objects, 
                                          "dapi_segment", 
                                          "dapi_internal"),
         nucleolus_internal = str_replace(objects, 
                                          "dapi_segment", 
                                          "nucleolus_internal"),
         nucleolus_external = str_replace(objects, 
                                          "dapi_segment", 
                                          "nucleolus_external"),
         nucleolus = str_replace(objects, 
                                 "dapi_segment", 
                                 "nucleolus_smooth"),
         tracer = str_replace(objects, 
                              "dapi_segment", 
                              "m6ATracer_smooth")) %>%
  arrange(cell, target, condition, replicate, slide)

# Print metadata
metadata %>% 
  print(n = Inf)

## # A tibble: 73 × 13
##    cell_statistics     sample  replicate cell  condition target slide objects   
##    <chr>               <chr>   <fct>     <fct> <chr>     <fct>  <chr> <chr>     
##  1 ../ts201102_E1349_… RPE_Da… r1        RPE   wt        Dam    002   ../ts2011…
##  2 ../ts201102_E1349_… RPE_Da… r1        RPE   wt        Dam    004   ../ts2011…
##  3 ../ts201102_E1349_… RPE_Da… r1        RPE   wt        Dam    006   ../ts2011…
##  4 ../ts201102_E1349_… RPE_Da… r1        RPE   wt        Dam    009   ../ts2011…
##  5 ../ts201102_E1349_… RPE_Da… r1        RPE   wt        Dam    011   ../ts2011…
##  6 ../ts201102_E1349_… RPE_Da… r1        RPE   wt        Dam    013   ../ts2011…
##  7 ../ts201102_E1349_… RPE_Da… r1        RPE   wt        Dam    015   ../ts2011…
##  8 ../ts201102_E1349_… RPE_Da… r1        RPE   wt        Dam    017   ../ts2011…
##  9 ../ts201102_E1349_… RPE_Da… r1        RPE   wt        Dam    019   ../ts2011…
## 10 ../ts210308_E1499_… RPE_Da… r2        RPE   wt        Dam    004   ../ts2103…
## 11 ../ts210308_E1499_… RPE_Da… r2        RPE   wt        Dam    007   ../ts2103…
## 12 ../ts210308_E1499_… RPE_Da… r2        RPE   wt        Dam    009   ../ts2103…
## 13 ../ts210308_E1499_… RPE_Da… r2        RPE   wt        Dam    012   ../ts2103…
## 14 ../ts210308_E1499_… RPE_Da… r2        RPE   wt        Dam    015   ../ts2103…
## 15 ../ts210308_E1499_… RPE_Da… r3        RPE   wt        Dam    002   ../ts2103…
## 16 ../ts210308_E1499_… RPE_Da… r3        RPE   wt        Dam    005   ../ts2103…
## 17 ../ts210308_E1499_… RPE_Da… r3        RPE   wt        Dam    007   ../ts2103…
## 18 ../ts210308_E1499_… RPE_Da… r3        RPE   wt        Dam    009   ../ts2103…
## 19 ../ts210308_E1499_… RPE_Da… r3        RPE   wt        Dam    011   ../ts2103…
## 20 ../ts201102_E1349_… RPE_Ki… r1        RPE   wt        Ki67   009   ../ts2011…
## 21 ../ts201102_E1349_… RPE_Ki… r1        RPE   wt        Ki67   011   ../ts2011…
## 22 ../ts201102_E1349_… RPE_Ki… r1        RPE   wt        Ki67   014   ../ts2011…
## 23 ../ts201102_E1349_… RPE_Ki… r1        RPE   wt        Ki67   016   ../ts2011…
## 24 ../ts201102_E1349_… RPE_Ki… r1        RPE   wt        Ki67   018   ../ts2011…
## 25 ../ts201102_E1349_… RPE_Ki… r1        RPE   wt        Ki67   021   ../ts2011…
## 26 ../ts201102_E1349_… RPE_Ki… r1        RPE   wt        Ki67   023   ../ts2011…
## 27 ../ts201102_E1349_… RPE_Ki… r1        RPE   wt        Ki67   025   ../ts2011…
## 28 ../ts201102_E1349_… RPE_Ki… r1        RPE   wt        Ki67   027   ../ts2011…
## 29 ../ts201102_E1349_… RPE_Ki… r1        RPE   wt        Ki67   029   ../ts2011…
## 30 ../ts201102_E1349_… RPE_Ki… r1        RPE   wt        Ki67   033   ../ts2011…
## 31 ../ts201102_E1349_… RPE_Ki… r1        RPE   wt        Ki67   035   ../ts2011…
## 32 ../ts210308_E1499_… RPE_Ki… r2        RPE   wt        Ki67   004   ../ts2103…
## 33 ../ts210308_E1499_… RPE_Ki… r2        RPE   wt        Ki67   006   ../ts2103…
## 34 ../ts210308_E1499_… RPE_Ki… r2        RPE   wt        Ki67   008   ../ts2103…
## 35 ../ts210308_E1499_… RPE_Ki… r2        RPE   wt        Ki67   010   ../ts2103…
## 36 ../ts210308_E1499_… RPE_Ki… r2        RPE   wt        Ki67   013   ../ts2103…
## 37 ../ts210308_E1499_… RPE_Ki… r2        RPE   wt        Ki67   015   ../ts2103…
## 38 ../ts210308_E1499_… RPE_Ki… r3        RPE   wt        Ki67   002   ../ts2103…
## 39 ../ts210308_E1499_… RPE_Ki… r3        RPE   wt        Ki67   004   ../ts2103…
## 40 ../ts210308_E1499_… RPE_Ki… r3        RPE   wt        Ki67   006   ../ts2103…
## 41 ../ts210308_E1499_… RPE_Ki… r3        RPE   wt        Ki67   008   ../ts2103…
## 42 ../ts210308_E1499_… RPE_Ki… r3        RPE   wt        Ki67   010   ../ts2103…
## 43 ../ts210308_E1499_… RPE_Ki… r3        RPE   wt        Ki67   012   ../ts2103…
## 44 ../ts210308_E1499_… RPE_Ki… r3        RPE   wt        Ki67   015   ../ts2103…
## 45 ../ts201008_E1304_… HCT116… r1        HCT1… DMSO      Dam    003   ../ts2010…
## 46 ../ts201008_E1304_… HCT116… r1        HCT1… DMSO      Dam    005   ../ts2010…
## 47 ../ts201008_E1304_… HCT116… r1        HCT1… DMSO      Dam    007   ../ts2010…
## 48 ../ts201008_E1304_… HCT116… r1        HCT1… DMSO      Dam    009   ../ts2010…
## 49 ../ts201008_E1304_… HCT116… r1        HCT1… wt        Dam    003   ../ts2010…
## 50 ../ts201008_E1304_… HCT116… r1        HCT1… wt        Dam    006   ../ts2010…
## 51 ../ts201008_E1304_… HCT116… r1        HCT1… wt        Dam    008   ../ts2010…
## 52 ../ts201008_E1304_… HCT116… r1        HCT1… wt        Dam    010   ../ts2010…
## 53 ../ts201008_E1304_… HCT116… r1        HCT1… wt        Dam    012   ../ts2010…
## 54 ../ts201008_E1304_… HCT116… r1        HCT1… DMSO      Ki67   003   ../ts2010…
## 55 ../ts201008_E1304_… HCT116… r1        HCT1… DMSO      Ki67   005   ../ts2010…
## 56 ../ts201008_E1304_… HCT116… r1        HCT1… DMSO      Ki67   007   ../ts2010…
## 57 ../ts201008_E1304_… HCT116… r1        HCT1… DMSO      Ki67   008   ../ts2010…
## 58 ../ts201008_E1304_… HCT116… r1        HCT1… DMSO      Ki67   009   ../ts2010…
## 59 ../ts201008_E1304_… HCT116… r1        HCT1… DMSO      Ki67   011   ../ts2010…
## 60 ../ts201008_E1304_… HCT116… r1        HCT1… DMSO      Ki67   012   ../ts2010…
## 61 ../ts201008_E1304_… K562_D… r1        K562  wt        Dam    005   ../ts2010…
## 62 ../ts201008_E1304_… K562_D… r1        K562  wt        Dam    007   ../ts2010…
## 63 ../ts201008_E1304_… K562_D… r1        K562  wt        Dam    009   ../ts2010…
## 64 ../ts201008_E1304_… K562_D… r1        K562  wt        Dam    011   ../ts2010…
## 65 ../ts201008_E1304_… K562_D… r1        K562  wt        Dam    013   ../ts2010…
## 66 ../ts201008_E1304_… K562_D… r1        K562  wt        Dam    015   ../ts2010…
## 67 ../ts201008_E1304_… K562_K… r1        K562  wt        Ki67   004   ../ts2010…
## 68 ../ts201008_E1304_… K562_K… r1        K562  wt        Ki67   006   ../ts2010…
## 69 ../ts201008_E1304_… K562_K… r1        K562  wt        Ki67   008   ../ts2010…
## 70 ../ts201008_E1304_… K562_K… r1        K562  wt        Ki67   010   ../ts2010…
## 71 ../ts201008_E1304_… K562_K… r1        K562  wt        Ki67   014   ../ts2010…
## 72 ../ts201008_E1304_… K562_K… r1        K562  wt        Ki67   016   ../ts2010…
## 73 ../ts201008_E1304_… K562_K… r1        K562  wt        Ki67   018   ../ts2010…
## # … with 5 more variables: cell_internal <chr>, nucleolus_internal <chr>,
## #   nucleolus_external <chr>, nucleolus <chr>, tracer <chr>

2) Load cell sizes and create a cell table

Load the cells and show the sizes. Possibly filter on those.

# Load all cells
LoadCells <- function(i) {
  tib <- read_csv(metadata$cell_statistics[i]) %>%
    dplyr::select(c(1:2)) %>%
    rename_all(~c("volume", "surface")) %>%
    add_column(object = 1:nrow(.),
               sample = metadata$sample[i],
               replicate = metadata$replicate[i],
               cell = metadata$cell[i],
               condition = metadata$condition[i],
               target = metadata$target[i],
               slide = metadata$slide[i]) %>%
    mutate(cell_id = paste(sample, object, sep = "_"))
  
  tib
  
}

tib_cells <- lapply(1:nrow(metadata), LoadCells) %>%
  purrr::reduce(bind_rows)

Plot and filter these cells.

# Set size limits
limits_size <- c(10, 125)

# Plot cell volume
tib_cells %>%
  ggplot(aes(x = target, y = volume)) +
  geom_hline(yintercept = 0, alpha = 0) +
  geom_hline(yintercept = limits_size, col = "black", linetype = "dashed") +
  geom_quasirandom(col = "darkgrey") +
  geom_boxplot(fill = NA, col = "black", outlier.shape = NA) +
  facet_grid(. ~ cell) +
  xlab("") +
  ylab("Volume slice (micron)") +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

# Filter cells for cutoff
tib_cells <- tib_cells %>%
  mutate(size_filter = volume >= limits_size[1] & 
           volume <= limits_size[2]) 

3) Load the pgm maps

Load the .pgm files for the various channels.

# Load all cells
LoadPGMMaps <- function(i, n_min = 5) {
  
  # Load objects
  objects <- ReadPGN(metadata$objects[i])
  nucleolus <- ReadPGN(metadata$nucleolus[i])
  nucleolus_internal <- ReadPGN(metadata$nucleolus_internal[i])
  nucleolus_external <- ReadPGN(metadata$nucleolus_external[i])
  tracer <- ReadPGN(metadata$tracer[i])
  
  tib <- tibble(object = c(objects),
                nucleolus_internal = c(nucleolus_internal),
                nucleolus_external = c(nucleolus_external),
                nucleolus = c(nucleolus),
                tracer = c(tracer)) %>%
    filter(object != 0) %>%
    rowwise() %>%
    mutate(nucleolus_internal = min(nucleolus_internal, 254),
           nucleolus_external = min(nucleolus_external, 254)) %>%
    ungroup() %>%
    group_by(object, nucleolus_internal, nucleolus_external) %>%
    summarise(n = n(),
              tracer = mean(tracer),
              nucleolus = mean(nucleolus)) %>%
    ungroup() %>%
    group_by(object) %>%
    mutate(tracer_enrichment = log2(tracer / weighted.mean(tracer, n)),
           nucleolus_enrichment = log2(nucleolus / weighted.mean(nucleolus, n))) %>%
    ungroup() %>%
    filter(n > n_min) %>%
    mutate(distance = nucleolus_external - nucleolus_internal) %>%
    arrange(object, distance) %>%
    add_column(sample = metadata$sample[i],
               replicate = metadata$replicate[i],
               cell = metadata$cell[i],
               condition = metadata$condition[i],
               target = metadata$target[i],
               slide = metadata$slide[i]) %>%
    mutate(cell_id = paste(sample, object, sep = "_"))
  
  tib
  
}

tib_intensities <- lapply(1:nrow(metadata), LoadPGMMaps) %>%
  purrr::reduce(bind_rows)

# Filter out cells that did not pass the cell size filter
tib_intensities <- tib_intensities %>%
  dplyr::filter(cell_id %in% (tib_cells %>% 
                                filter(size_filter == T) %>%
                                pull(cell_id)))

Plots.

# Filter very bad cells
# - only positive / negative scores
# - extreme values
remove_samples <- tib_intensities %>% 
  gather(key, value, nucleolus_enrichment, tracer_enrichment) %>%
  mutate(key = str_remove(key, "_enrichment")) %>%
  filter(distance > -11 & distance <= 11) %>% 
  group_by(key, target, cell, cell_id) %>% 
  dplyr::summarise(q1 = quantile(value, 0.1), 
                   q5 = quantile(value, 0.5), 
                   q9 = quantile(value, 0.9)) %>%
  filter(q1 > 0.5 | q9 < 0 | q1 < -2.5 | q9 > 3) %>%
  pull(cell_id)

## `summarise()` has grouped output by 'key', 'target', 'cell'. You can override using the `.groups` argument.

# Plot intensity per cell
tib_gather <- tib_intensities %>%
  filter(! cell_id %in% remove_samples) %>%  
  gather(key, value, nucleolus, tracer) %>%
  filter(distance > -14 & distance <= 14) %>% 
  mutate(key = factor(key, levels = c("tracer", "nucleolus")))

tib_gather %>%  
  ggplot(aes(x = distance, y = value, col = key, 
             group = interaction(key, cell_id, condition))) +
  geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
  geom_vline(xintercept = 0, col = "black", linetype = "dashed") +
  geom_line() +
  facet_grid(target ~ cell, scales = "free_y") +
  xlab("Distance to nucleolus (pixel)") +
  ylab("Intensity (A.U.)") +
  coord_cartesian(xlim = c(-10, 10)) +
  scale_color_grey() +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

# Plot intensity per cell - for enrichments
tib_gather <- tib_intensities %>%
  filter(! cell_id %in% remove_samples) %>%  
  gather(key, value, nucleolus_enrichment, tracer_enrichment) %>%
  filter(distance > -14 & distance <= 14) %>% 
  mutate(key = str_remove(key, "_enrichment")) %>%
  mutate(key = factor(key, levels = c("tracer", "nucleolus")))

tib_gather %>%  
  ggplot(aes(x = distance, y = value, col = key, 
             fill = target, group = interaction(key, cell_id, condition))) +
  geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
  geom_vline(xintercept = 0, col = "black", linetype = "dashed") +
  geom_line(alpha = 0.1) +
  stat_summary(fun = mean, aes(group = key), 
               geom = "line", size = 2) +
  facet_grid(cell ~ target, scales = "free_y") +
  xlab("Distance to nucleolus (pixel)") +
  ylab("Enrichment over mean (log2)") +
  coord_cartesian(xlim = c(-10, 10),
                  ylim = c(-2.1, 2.3)) +
  scale_color_manual(values = c("grey30", "red")) +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

tib_gather %>%  
  ggplot(aes(x = distance, y = value, col = target, 
             fill = target, 
             group = interaction(target, cell_id, condition))) +
  geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
  geom_vline(xintercept = 0, col = "black", linetype = "dashed") +
  geom_line(alpha = 0.1) +
  stat_summary(fun = mean, aes(group = target), 
               geom = "line", size = 2) +
  # stat_summary(fun.data = mean_se, geom = "ribbon", aes(group = target),
  #              fun.args = list(mult = 1.96), col = NA, alpha = 0.2) +
  facet_grid(cell ~ key, scales = "free_y") +
  xlab("Distance to nucleolus (pixel)") +
  ylab("Enrichment over mean (log2)") +
  coord_cartesian(xlim = c(-10, 10),
                  ylim = c(-2.1, 2.3)) +
  scale_color_manual(values = c("grey30", "red")) +
  scale_fill_manual(values = c("grey30", "red")) +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

# Means only
# 1) Combined
tib_gather %>%  
  ggplot(aes(x = distance, y = value, col = target, 
             fill = target, group = target)) +
  geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
  geom_vline(xintercept = 0, col = "black", linetype = "dashed") +
  stat_summary(fun = mean, geom = "line", size = 1) +
  # stat_summary(fun.data = mean_se, geom = "errorbar",
  #              col = "darkgrey", linetype = "solid") +
  stat_summary(fun.data = mean_se, geom = "ribbon",
               col = NA, alpha = 0.2) +
  facet_grid(cell ~ key) +
  xlab("Distance to nucleolus (pixel)") +
  ylab("Enrichment over mean (log2)") +
  coord_cartesian(xlim = c(-10, 10),
                  ylim = c(-1.5, 1.9)) +
  #scale_color_grey() +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

# 2) Per replicate
tib_gather %>%  
  ggplot(aes(x = distance, y = value, col = target, linetype = interaction(replicate, condition),
             fill = target, group = interaction(target, replicate, condition))) +
  geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
  geom_vline(xintercept = 0, col = "black", linetype = "dashed") +
  stat_summary(fun = mean, geom = "line", size = 1) +
  # stat_summary(fun.data = mean_se, geom = "errorbar",
  #              col = "darkgrey", linetype = "solid") +
  stat_summary(fun.data = mean_se, geom = "ribbon",
               col = NA, alpha = 0.2) +
  facet_grid(cell ~ key) +
  xlab("Distance to nucleolus (pixel)") +
  ylab("Enrichment over mean (log2)") +
  coord_cartesian(xlim = c(-10, 10),
                  ylim = c(-1.5, 1.9)) +
  #scale_color_grey() +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

4) Distance maps for cells instead of nucleoli

As requested by reviewer #2, include quantification from periphery of nuclei

# Load all cells
LoadPGMMaps <- function(i, n_min = 5) {
  
  # Load objects
  objects <- ReadPGN(metadata$objects[i])
  cell_internal <- ReadPGN(metadata$cell_internal[i])
  nucleolus <- ReadPGN(metadata$nucleolus[i])
  tracer <- ReadPGN(metadata$tracer[i])
  
  tib <- tibble(object = c(objects),
                cell_internal = c(cell_internal),
                nucleolus = c(nucleolus),
                tracer = c(tracer)) %>%
    filter(object != 0) %>%
    rowwise() %>%
    mutate(cell_internal = min(cell_internal, 254)) %>%
    ungroup() %>%
    group_by(object, cell_internal) %>%
    summarise(n = n(),
              tracer = mean(tracer),
              nucleolus = mean(nucleolus)) %>%
    ungroup() %>%
    group_by(object) %>%
    mutate(tracer_enrichment = log2(tracer / weighted.mean(tracer, n)),
           nucleolus_enrichment = log2(nucleolus / weighted.mean(nucleolus, n))) %>%
    ungroup() %>%
    filter(n > n_min) %>%
    mutate(distance = 254 - cell_internal) %>%
    arrange(object, distance) %>%
    add_column(sample = metadata$sample[i],
               replicate = metadata$replicate[i],
               cell = metadata$cell[i],
               condition = metadata$condition[i],
               target = metadata$target[i],
               slide = metadata$slide[i]) %>%
    mutate(cell_id = paste(sample, object, sep = "_"))
  
  tib
  
}

tib_intensities <- lapply(1:nrow(metadata), LoadPGMMaps) %>%
  purrr::reduce(bind_rows)

# Filter out cells that did not pass the cell size filter
tib_intensities <- tib_intensities %>%
  dplyr::filter(cell_id %in% (tib_cells %>% 
                                filter(size_filter == T) %>%
                                pull(cell_id)))

Plots.

# Filter same bad cells as before


# Plot intensity per cell
tib_gather <- tib_intensities %>%
  filter(! cell_id %in% remove_samples) %>%  
  gather(key, value, nucleolus, tracer) %>%
  filter(distance <= 14) %>% 
  mutate(key = factor(key, levels = c("tracer", "nucleolus")))

tib_gather %>%  
  ggplot(aes(x = distance, y = value, col = key, 
             group = interaction(key, cell_id, condition))) +
  geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
  geom_vline(xintercept = 0, col = "black", linetype = "dashed") +
  geom_line() +
  facet_grid(target ~ cell, scales = "free_y") +
  xlab("Distance from periphery (pixel)") +
  ylab("Intensity (A.U.)") +
  coord_cartesian(xlim = c(0, 10)) +
  scale_color_grey() +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

# Plot intensity per cell - for enrichments
tib_gather <- tib_intensities %>%
  filter(! cell_id %in% remove_samples) %>%  
  gather(key, value, nucleolus_enrichment, tracer_enrichment) %>%
  filter(distance > -14 & distance <= 14) %>% 
  mutate(key = str_remove(key, "_enrichment")) %>%
  mutate(key = factor(key, levels = c("tracer", "nucleolus")))

tib_gather %>%  
  ggplot(aes(x = distance, y = value, col = key, 
             fill = target, group = interaction(key, cell_id, condition))) +
  geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
  geom_vline(xintercept = 0, col = "black", linetype = "dashed") +
  geom_line(alpha = 0.1) +
  stat_summary(fun = mean, aes(group = key), 
               geom = "line", size = 2) +
  facet_grid(cell ~ target, scales = "free_y") +
  xlab("Distance from periphery (pixel)") +
  ylab("Enrichment over mean (log2)") +
  coord_cartesian(xlim = c(0, 10),
                  ylim = c(-3, 0.5)) +
  scale_color_manual(values = c("grey30", "red")) +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

tib_gather %>%  
  ggplot(aes(x = distance, y = value, col = target, 
             fill = target, 
             group = interaction(target, cell_id, condition))) +
  geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
  geom_vline(xintercept = 0, col = "black", linetype = "dashed") +
  geom_line(alpha = 0.1) +
  stat_summary(fun = mean, aes(group = target), 
               geom = "line", size = 2) +
  # stat_summary(fun.data = mean_se, geom = "ribbon", aes(group = target),
  #              fun.args = list(mult = 1.96), col = NA, alpha = 0.2) +
  facet_grid(cell ~ key, scales = "free_y") +
  xlab("Distance from periphery (pixel)") +
  ylab("Enrichment over mean (log2)") +
  coord_cartesian(xlim = c(0, 10),
                  ylim = c(-3, 0.5)) +
  scale_color_manual(values = c("grey30", "red")) +
  scale_fill_manual(values = c("grey30", "red")) +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

# Means only
# 1) Combined
tib_gather %>%  
  ggplot(aes(x = distance, y = value, col = target, 
             fill = target, group = target)) +
  geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
  geom_vline(xintercept = 0, col = "black", linetype = "dashed") +
  stat_summary(fun = mean, geom = "line", size = 1) +
  # stat_summary(fun.data = mean_se, geom = "errorbar",
  #              col = "darkgrey", linetype = "solid") +
  stat_summary(fun.data = mean_se, geom = "ribbon",
               col = NA, alpha = 0.2) +
  facet_grid(cell ~ key) +
  xlab("Distance from periphery (pixel)") +
  ylab("Enrichment over mean (log2)") +
  coord_cartesian(xlim = c(0, 10),
                  ylim = c(-2, 0.2)) +
  #scale_color_grey() +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

# 2) Per replicate
tib_gather %>%  
  ggplot(aes(x = distance, y = value, col = target, linetype = interaction(replicate, condition),
             fill = target, group = interaction(target, replicate, condition))) +
  geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
  geom_vline(xintercept = 0, col = "black", linetype = "dashed") +
  stat_summary(fun = mean, geom = "line", size = 1) +
  # stat_summary(fun.data = mean_se, geom = "errorbar",
  #              col = "darkgrey", linetype = "solid") +
  stat_summary(fun.data = mean_se, geom = "ribbon",
               col = NA, alpha = 0.2) +
  facet_grid(cell ~ key) +
  xlab("Distance from periphery (pixel)") +
  ylab("Enrichment over mean (log2)") +
  coord_cartesian(xlim = c(0, 10),
                  ylim = c(-2, 0.2)) +
  #scale_color_grey() +
  theme_bw() +
  theme(aspect.ratio = 1,
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

Conclusion

Several observations:

• Ki67 m6A-Tracer is enriched near Ki67 protein. That’s good.
• This enrichment is almost equal between HCT116 and K562 cells. RPE cells are different, probably due to a difference in slide quality. At least, the nucleolus difference is also less pronounced at nucleolar borders.

Overall, this is a positive result.

SessionInfo

sessionInfo()

## R version 4.1.2 (2021-11-01)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS Monterey 12.0.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1-arm64/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] knitr_1.36       pixmap_0.4-12    ggbeeswarm_0.6.0 forcats_0.5.1   
##  [5] stringr_1.4.0    dplyr_1.0.7      purrr_0.3.4      readr_2.1.0     
##  [9] tidyr_1.1.4      tibble_3.1.6     tidyverse_1.3.1  ggplot2_3.3.5   
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.7       lubridate_1.8.0  assertthat_0.2.1 digest_0.6.28   
##  [5] utf8_1.2.2       R6_2.5.1         cellranger_1.1.0 backports_1.4.0 
##  [9] reprex_2.0.1     evaluate_0.14    highr_0.9        httr_1.4.2      
## [13] pillar_1.6.4     rlang_0.4.12     readxl_1.3.1     rstudioapi_0.13 
## [17] jquerylib_0.1.4  rmarkdown_2.11   labeling_0.4.2   bit_4.0.4       
## [21] munsell_0.5.0    broom_0.7.10     compiler_4.1.2   vipor_0.4.5     
## [25] modelr_0.1.8     xfun_0.28        pkgconfig_2.0.3  htmltools_0.5.2 
## [29] tidyselect_1.1.1 codetools_0.2-18 fansi_0.5.0      crayon_1.4.2    
## [33] tzdb_0.2.0       dbplyr_2.1.1     withr_2.4.2      grid_4.1.2      
## [37] jsonlite_1.7.2   gtable_0.3.0     lifecycle_1.0.1  DBI_1.1.1       
## [41] magrittr_2.0.1   scales_1.1.1     vroom_1.5.6      cli_3.1.0       
## [45] stringi_1.7.5    farver_2.1.0     fs_1.5.0         xml2_1.3.2      
## [49] ellipsis_0.3.2   generics_0.1.1   vctrs_0.3.8      tools_4.1.2     
## [53] bit64_4.0.5      glue_1.5.0       beeswarm_0.4.0   hms_1.1.1       
## [57] parallel_4.1.2   fastmap_1.1.0    yaml_2.2.1       colorspace_2.0-2
## [61] rvest_1.0.2      haven_2.4.3